C======================================================================      
C          This subroutine sets up the r-mesh.
C======================================================================      
      subroutine  gridr1 (hmin, emax, npwave, hmax, rmax, rdbl, 
     >   gridr, nr, idouble, ndouble)
      dimension  gridr(nr, 3), idouble(0:ndouble)
C
C  INPUT:
C  -----
C   hmin - the first step of the mesh.
C   Emax - the biggest energy what could be calculated by this mesh.
C   npwave - number of points per wave with Emax, more then 3.
C            Normal Npwave is 6.
C   rmax - the largest "r" in the mesh.
C   rdbl - the point of the last doubling.
C
C  OUTPUT:
C  ------
C   hmax - the last step of the mesh.
C   nr - total number of points in the mesh.
C   gridr(i,1) - r(i)      points. 
C   gridr(i,2) - hr(i)     steps. 
C   gridr(i,3) - wr(i)     weigths. 
C
C  Idouble stores the I index for GRIDR where DR doubles, with the first
C  point set to 1 and the last to NR Used in the numerov integrations.
C
C--------------------------------------------------------------------------
C
C     Do some check 
C
      if (npwave .lt. 4) then
         print*, 'Npwave is too little'
         stop    'Npwave is too little'
      end if 
      hmax = 3.14 / float(npwave) / sqrt(emax * 2.0)
      if ((hmin .le. 0.0) .or. (hmin .gt. hmax)) then
         print*, 'HMIN incompatible with HMAX:', hmin, hmax
         stop    'HMIN incompatible with HMAX'
      end if 
C      
C     2**ndouble * hmin = hmax, therefore
C
      ntmp = int(log(hmax/hmin)/log(2.0))
      if (rdbl .le. 0.0) then
         ntmp = 0
         hmax = hmin
      end if 
      if (ntmp .gt. ndouble) then
         print*, 'Increase NDOUBLE in call to GRIDS to at least:', ntmp
         stop    'Increase NDOUBLE in call to GRIDS'
      end if 
      ndouble = ntmp
      dr = hmin
C
C     let NPDBL be the number of points with the same dr per interval
C     RDBL = NPDBL * hmin * (2**NDOUBLE-1)
C
      if (ndouble .gt. 0) then
         npdbl = rdbl / (hmin * float(2**ndouble-1))
      else
         npdbl = 0
      end if 
C     
C     Make sure NPDBL is even
C
      npdbl=(npdbl / 2) * 2 
      if (npdbl .lt. 4 .and. npdbl .ne. 0) then
         print*,'Need to have at least 4 equally spaced points in GRIDR'
         stop   'Need to have at least 4 equally spaced points in GRIDR'
      end if 
      idouble(1) = 1
      do nd = 2, ndouble + 1
         idouble(nd) = npdbl * (nd - 1)
      end do 
C
C     Set up r-mesh 
C
      r = 0.0
      i = 0
      do nd = 1, ndouble
         do ni = 1, npdbl
            i = i + 1
            gridr(i, 1) = r + float(ni) * dr
            gridr(i, 2) = dr
C
C           Simpson's rule weights
C
            gridr(i, 3) = float(mod(i, 2) * 2 + 2) * dr / 3.0
         end do
         r = gridr(i, 1)
         gridr(i, 3) = dr
         dr = dr * 2.0
      end do
c$$$      ni = nint((rmax-r)/hmax)
      if (dr .gt. hmax)    dr = dr / 2.0
      hmax = dr 
      ni = nint((rmax - r) / dr)
      do ii = 1, ni
         i = i + 1
         gridr(i, 1) = r + float(ii) * dr
         gridr(i, 2) = dr
C
C        Simpson's rule weights
C
         gridr(i, 3) = float(mod(i, 2) * 2 + 2) * dr / 3.0
      end do
      if (i.gt.nr) then
         print*, 'Increase NR in call to GRIDS to at least:', i
         stop    'Increase NR in call to GRIDS'
      end if 
      nr = i
      idouble(ndouble + 2) = i
      ntmp = nr 
      if (ntmp .gt. nr) then
         print*, 'Increase NR in call to GRIDS to at least:', ntmp
         stop    'Increase NR in call to GRIDS'
      end if 
      return
      end
c
      subroutine weights_c (T, WTS, NT, tmp1, tmp2, tmp3)
      integer NT, i
      real*8 T(NT), WTS(NT), tmp1, tmp2, tmp3
      do i = 1, NT
         T(i)   = tmp1 * T(i) - tmp2
         WTS(i) = tmp3 * WTS(i)
      end do 
      return
      end
c
